Elastic constants of /3-eucryptite: A density functional theory study 
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The five independent elastic constants of hexagonal /3-eucryptite have been determined using den- 
sity functional theory (DFT) total energy calculations. The calculated values agree well, to within 
15%, with the experimental data. Using the calculated elastic constants, the linear compressibility 
of /3-eucryptite parallel to the c-axis, \c, and perpendicular to it, \ a , have been evaluated. These 
values are in close agreement to those obtained from experimentally known elastic constants, but 
are in contradiction to the direct measurements based on a three-terminal technique. The calculated 
compressibility parallel to the c-axis was found to positive as opposed to the negative value obtained 
by direct measurements. We demonstrate that Xc must be positive and discussed the implications 
of a positive Xc in the context of explaining the negative bulk thermal expansion of /3-eucryptite. 



I. INTRODUCTION 

Glass ceramics in Li 2 0-Al 2 03-Si02 (LAS) systems 
have attracted a lot of attention over the last several 
decades due to their low or even negative coefficient of 
thermal expansion (CTE), as well as due to their chem- 
ical and thermal stability^- 3 - This class of materials has 
been extensively commercialized owing to their exotic 
physical properties which makes them suitable for indus- 
trial applications (e.g., heat exchangers) which require 
dimensional stability and thermal shock resistances^ 
They are also used in very specific applications like tele- 
scope mirror blanks, high precision optical devices and 
ring laser gyroscope^! 3 - The hexagonal /3-eucryptite is 
a prominent member of this class of materials. It has 
a highly anisotropic CTE 3 (i.e., a a = 7.26 xlO -6 per- 
pendicular to the c axis, a c = -16.35 xlCT 6 parallel to 
the c axis) which leads to a slightly negative crystallo- 
graphic average (bulk) CTE. /3-eucryptite undergoes a 
reversible order-disorder structural transition at ~ 755 
Ki£ It exhibits one dimensional superionic conductivity 
of Li + ions along the c-axis which makes it a suitable 
electrolyte in Li based batteries^ Most of these unusual 
properties of /3-eucryptite are, in part, related to its crys- 
tal structure. 

Figure Q] illustrates a unit cell of /3-eucryptite below 
the order-disorder transition containing 84 atoms with 
12 unit formulae of LiAlSi04. A single crystal of ordered 
/3-eucryptite, as shown by Figure^has a primitive hexag- 
onal structure belonging to the P6422 space group.— This 
structure is a derivative of the /3-quartz configuration, 
with half the Si 4+ ions replaced by Al 3+ while the charge 
imbalance is compensated by the channels of Li + ions 
parallel to the c-axis.— Several researchers 3 -^^— have 
demonstrated through structural refinements that the 
structure is composed of interconnected helices of Si04 4 ~ 
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and AIO4 tetrahedra with alternation of layers contain- 
ing Si and Al atoms respectively, leading to a doubling 
of the c-axis of /3-quartz. 

As mentioned, the slightly negative crystallographic 



FIG. 1: (Color online) Crystal structure of ordered ,9- 
eucryptite containing 12 formula units of LiAlSi04 (84 
atoms) . 



average CTE of /3-eucryptite is due to the anisotropy 
of the linear expansion where a temperature increase in- 
volves a contraction along the c-axis which overcompen- 
sates the concomitant expansion in the plane perpen- 
dicular to the c-axis. Several theories have been pro- 
posed to explain this unusual thermal behaviori 3 -^ - — 
Hortal et alJ^- employed a three-terminal technique 1 ^ 
to measure the linear compressibility \ of /3-eucryptite 
along the a and c axes, and the reported values of Xa — 
(22.4±6.0) x 1CT 3 GPa" 1 and X c = (-1.13±1.0) x 1(T 3 
GPa^ 1 . These measured values of \ supported the ex- 
planation given by Gillery and Bush 1 — that the negative 
bulk thermal expansion is an elastic effect associated with 
the interconnected helices of Si and Al tetrahedra. 
It is well known that linear compressibility along any 
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direction in a crystal can be calculated from the elements 
of the stiffness matrix CV,-. 19 Linear compressibilities 
can also be calculated from experimentally determined 20 
stiffness constants; we carried out this calculation and the 
corresponding error analysis, and have found that Xa — 
(2.57±0.02) x 1(T 3 GPa 11 and X c = (4.60±0.05) x 1(T 3 
GPa -1 at a temperature of 293K. It should be noted 
that Xc calculated using the experimentally known Cjj 
is positive, in contrast to the negative value of (—1.13 ± 
1.0) x 10~ 3 GPa determined by direct measurements.— 
While the sign of the direct measurement of Xc remains 
in doubt due to experimental uncertainty, the near-zero 
value is also very different from that obtained in calcula- 
tions. Since the sign of Xc is linked with the explanation 
of negative CTE of /3-eucryptite, it is necessary to ad- 
dress and possibly resolve the contradiction surrounding 
the sign of Xc- 

In this paper, we compute the elastic stiffness con- 
stants of ordered /3-eucryptite containing 84 atoms per 
unit cell in the framework of density functional theory 
(DFT) . We then use the elastic constants to evaluate the 
linear compressibilities Xa & n d Xc to clarify the sign of 
linear compressibility parallel to the c-axis. Since the 
density functional theory (DFT) calculations offer an in- 
dependent method of determining the linear compress- 
ibility, the present study can resolve the discrepancy dis- 
cussed above. After ascertaining the sign of Xc, we dis- 
cuss its implications on the explanation of negative crys- 
tallographic average CTE of /3-eucryptite. We demon- 
strate that the negative CTE of /3-eucryptite must arise 
from a combination of several interconnected phenomena 
as suggested by Xu et al£ and is related to a negative 
Griineissen function along the c-axis, rather than to the 
elastic effect proposed by Gillery and Bush. 15 

The paper is organized as follows: Sec II describes the 
methodology we adopted to calculate the elastic stiffness 
constants and the details of the DFT calculations; Sec 
III describes the results obtained in the present study 
which are discussed in Sec IV in the context of resolving 
the discrepancy and explaining the negative CTE of /3- 
eucryptite; Sec V summarizes the results and describes 
our main conclusions. 



II. METHODOLOGY 
A. Calculation of elastic constants 

In general, a crystal deforms in a homogeneous lin- 
ear elastic manner when subjected to sufficiently small 
strains (i,j = 1,2,3). The components Cyy of the 
adiabatic stiffness matrix are the derivatives of elastic 
energy density with respect to the strain components: 21 



where E is the elastic energy stored in a domain of volume 
V of the crystal subjected to homogeneous deformations, 
and Vq is the volume of the unstrained crystal. 

In this section, we briefly describe the technique 22 
we employed to calculate the elastic constants of (3- 
eucryptite. The lattice of hexagonal /3-eucryptite is 
spanned by three primitive Bravais lattice vectors which 
can be written in a matrix form as: 
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where each row is a lattice vector, and a, c are the two 
lattice parameters that characterize the hexagonal struc- 
ture. The vectors of the deformed lattice (R') can be 
obtained by multiplying R with a distortion matrix D: 



R' = RD, 



(3) 



where D is defined in terms of the components of strain 
tensor as: 
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(4) 



The elastic energy E of a crystal subjected to a general 
elastic strain (e^) can be expressed by means of a Taylor 
expansion in the distortion parameters truncated at the 
second order of strain, 21 
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where Eq is the energy of a crystal volume Vq at equi- 
librium, <jjj are the elements of the stress tensor, and 5ij 
is the Kronecker symbol. Since the distortion matrix is 
symmetric, it is convenient to express Eq. ([5]) using the 
Voigt notation (11 = 1, 22 = 2, 33 = 3, 23 = 4, 31 = 5, 
and 12 = 6): 



E(V,e) = E + V 



(6) 



where r\\ = 1 if i = 1, 2, or 3 and rji — 2 if i = 4, 5, or 6. 

Due to the specific symmetry of the hexagonal lattice, 
there are only five independent elastic constants, 19 which 
in Voigt notation are Cn, C\%, C13, C33 and C44. These 
constants can be determined from specific distortion ma- 
trices for the hexagonal structures^ Table U lists the 
distortion matrices used in the present study; for these 
matrices Eq. (jBJ) takes the simple form 



E(V,5) =E + V (A 1 5 + A 2 S 2 ), 



(7) 
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where A\ is related to stress components Uij , and Ai is a 
linear combination of the elastic constants Gy- The re- 
lationships between the second-order coefficient A2 and 
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TABLE I: The distortion matrices and elastic constants for a 
hexagonal lattice. 

Distortion matrix Second-order coefficient 

A 2 in Eq. © 



Cn + C12 




1 + 6 
1 + 8 
1 + 5 



Cn — Cl2 



Cn + Cl2 



+2C 13 + ifa 




C 33 /2 



2C 4 



the independent elastic constants for different strains are 
also given in Table Q] For each of the five different defor- 
mations listed in Table HI the total energy of the crystal 
was calculated for values of parameter S varying from 
—5% to +5%. The zeroth, first, and second order coeffi- 
cients in Eq. ([7]) were extracted by means of polynomial 
fits of the total energy versus 6 data. Using the relation- 
ships in the left column of Table []] the elastic constants 
were extracted from the coefficients A2 of the distortions 
considered. For all the cases, we found that the contribu- 
tion of order 3 and higher terms to the energy in Eq. (|5|) 
was negligible, which confirms that the strains used are 
within the linear elastic limit. 



B. Details of density functional calculations 

The total energy calculations for /3-eucryptite super- 
cells were performed within the framework of density 
functional theory (DFT) using the generalized gradi- 
ent approximation (GGA) and the projector augmented 
wave (PAW) potentials 23 implemented in the ab-initio 
simulation package VASP. 2 -^ We used the Perdew- 
Wang exchange-correlation function.— The plane wave 
energy cutoff was set to 500 eV, and the augmentation 
charge cutoff to 605 eV. The computational cell consisted 



of one primitive cell of hexagonal /3-eucryptite having 84 
atoms, i.e. 12 formula units (f.u.) of LiAlSiO.4. The 
Brillouin zone was sampled with a 3x3x3 T-centered 
Monkhorst-Pack grid, which resulted in 14 irreducible k- 
points. The atomic coordinates were optimized using the 
conjugate gradient algorithm until the force components 
on any atom were smaller than 0.01 eV/A. 



III. RESULTS 

To determine the equilibrium lattice parameters a and 
c of /3-eucryptite, the ab-initio total energy calculations 
were performed for different values of the supercell vol- 
ume V and the c/a ratio. The volume of the supercell was 
varied from —15% to +15% of the experimental value&iS 
by changing the parameter a while keeping the ratio c/a 
fixed. The procedure for finding the lattice constants a 
and c can be followed on the schematics in Fig. [Ha). At 
any given volume V [which is a const ant- volume plane 
in Fig. EJa)], the supercell was relaxed and the total en- 
ergy E was computed for different c/a values ranging 
from 1.00 to 1.12; the E vs. c/a data was then fitted 
with a fourth-order polynomial to determine the mini- 
mum energy for the given volume V. These minimum 
energy values found for different volumes V are plotted 
in Fig. HJb) and show an excellent fit to the Murnaghan 
equation of stated 



E(V) = E Q 
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where Eq is the energy corresponding to the equilibrium 
volume Vo, B is the bulk modulus at zero pressure, and 
B' = (OB /dP) T=Q is the pressure derivative of the bulk 
modulus at K. The parameters obtained from the fitting 
of E versus V data against the Murnaghan equation of 
state are listed in Table ITfl 

The bulk modulus obtained from the fit, as shown by 
Table UH is in excellent agreement with the experimen- 
tally determined valued which indicates a good agree- 
ment of the calculated E versus V data with Eq. (JSJ). 
At the equilibrium volume Vo obtained from the Eq. (JSJ , 
the energy of the supercell was calculated for different 
c/a values and plotted in Fig. [2jc). The E vs. c/a data 
at constant volume Vq was fitted against a fourth-order 
polynomial to determine the optimum c/a ratio at Vo- 
The optimized lattice constants a and c are subsequently 
extracted from the optimum c/a ratio and the equilib- 
rium volume Vo- The calculated values a and c (see Ta- 
ble [n]) are in close agreement to the experimental ones 
a exp and c exp , i.e., a = l.Qla exp and c = 1.015c exp , which 
lead to V = 1.03V exp . 

Using the calculated lattice parameters, we determined 
the five independent elastic constants of /3-eucryptite at 
K by employing the technique outlined in Section II. Fig- 
ure [3] shows the energy as a function of the deformation 
parameter 6 for the different types of strain listed in Ta- 
ble U along with the corresponding fit polynomials. The 
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FIG. 2: Determination of lattice constants a and c from total 
energy GGA calculations, (a) Schematic representation of the 
procedure to obtain the structural parameters of /3-eucryptite 
at equilibrium, (b) Energy vs. volume curve for the supercell 
shown in Fig. Q] The black squares are the DFT calculated 
points, while the solid line represents the fit to the Murnaghan 
equation of state, (c) Energy as a function of the ratio c/a at 
the equilibrium volume Vo given by the Murnaghan fit. The 
squares are the DFT calculated values, and the solid line is a 
fourth-order polynomial fit. 



elastic constants were evaluated from the second order 
coefficients of the fit polynomials through their rela- 
tionships with the stiffness constants listed in Table HI 

The elastic constants have been measured experimen- 
tally by Haussuhl et al. using an ultrasonic technique 
at ambient temperature, 293K* 20 In order to compare at 
the same temperature the values of Cij computed in the 



TABLE II: Calculated lattice parameters, equilibrium vol- 
ume Vo, bulk modulus B, and its pressure derivative (£>' = 
(dB /dP) T=0 ) of /3-eucryptite. The experimental values are 
also provided wherever available. 



Technique q(A) c(A) V (A 3 /uf) B(GPa) B' 

GGA 10.594 11.388 92.25 102.27 -1.05 

Exp 10.497° 11.200" 89.06" 109.9 6 



a Ref.[Tcl 293K. 
fc Prom Cij in Ref. 
extrapolated to K. 



TABLE III: Comparison of the calculated stiffness constants 
Cij of /3-eucryptite with the experimental data from Ref. 
extrapolated to K using the thermoelastic constants Tij = 
dlog dj/dT. The uncertainty in any of the experimental val- 
ues (Exp) is smaller than 2.5 GPa. 



Cn Cl2 C13 C33 C44 



GGA: (GPa) 165.64 70.98 78.59 132.83 58.68 
Exp: (GPa) 176.3 68.5 89.8 139.9 61.2 
Tij {W~ 3 /K) -0-14 0.13 -0.27 -0.42 -0.24 



present study (GGA) with the experiments, we have ex- 
trapolated the measured values of Cy to K by using 
the thermoelastic constants = d log Cy /dT*& The 
calculated clastic constants that we obtained are within 
^15% of the experimental values extrapolated to K (see 
Table ITTT|) . Furthermore, we have also found that the ex- 
trapolation of the GGA elastic constants to 273 K is also 
consistent with the experimental data at 273 K. 

The linear compressibilities Xa and \ c (along the a and 
c axes) for a transversely isotropic material are related 
to the elastic constants dj through: 19 

C33 — C13 

Xa 0/^2 1 /~< 

1^11^33 — z< ~-13 + L , 12<^33 

C\\ + C12 — 2C13 / q x 

L'llL'SS — ^13 + ^12^33 

With the calculated elastic constants (Table IIIip in 
Eqs. ©, we determined the linear compressibilities of 
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FIG. 3: Total energy of ordered /3-eucryptite as a function of 
deformation parameter 8 for the five different distortions (a) 
- (e) in the same order as listed in Table [I] 



/3-eucryptite at K. However, the corresponding experi- 
mental data- obtained by direct measurements using a 
three-terminal method— have been reported at 273 K. 
In order to make comparisons of compressibility values 
at the same temperature we have extrapolated the CV, 
values from GGA calculations to 273 K using the ther- 
moelastic constants from Rcf. |2y. Similarly, we eval- 
uated the experimentally determined elastic constants at 
273 K. The extrapolated values of the elastic constants 
were then used in Eq. (J9j> to determine Xa and Xc, which 
were compared with the direct measurements fTablc llVl) . 
We should note that the uncertainty in the GGA values 



TABLE IV: Comparison of calculated values of linear com- 
pressibility of ordered /3-eucryptite \a an d Xc with experi- 
mental data. The GGA elastic constants and the experimen- 
tal values-^ have been extrapolated to 273 K to calculate \a 
and Xc The calculated values are in good agreement with 
each other. However, they are in contradiction to the direct 
measurements^ at 273 K. 



Parameter GGA Ref. 20 Ref. 17 
(10~ 3 GPa' 1 ) (10~ 3 GPa' 1 ) (10~ 3 GPa" 1 ) 

Xa 2.67 ± 0.06 2.58 ± 0.02 22.4 ± 6.0 

Xc 5.20 ± 0.15 4.52± 0.05 -1.13 ± 1.0 



for compressibility comes solely from propagating the un- 
certainties in the thermoelastic constants^ 



IV. DISCUSSION 

Our calculated linear compressibility values of or- 
dered /3-eucryptite, extrapolated to 273 K, agree well 
with those derived using experimentally determined 
stiffness constants^ These two sets of compressibility 
values, however, are in contradiction with the direct 
measurements^ reported by Hortal et alJ^- We note that 
the calculated values of Xc are positive, while the value 
reported from direct measurements is negative (refer to 
Table IIV[) . Furthermore, the measured value of Xa is 
about one order of magnitude larger than that calculated 
in the present study. 

We now focus on the implications of the calculated 
compressibility values on the thermal behavior of /3- 
eucryptite, in particular on the coefficient of thermal ex- 
pansion. An early study by Munn^S addresses the ef- 
fect of anisotropy of elastic properties on the thermal 
expansion in the quasi-harmonic approximation, where 
the the vibrations are taken to be harmonic but with 
deformation-dependent frequencies. The bulk thermal 
expansion coefficient a of a hexagonal crystal can be writ- 
ten in terms of the thermal expansion coefficients along 
the c-axis (a c ) and a-axis (a a ) as^ 



2a a 



3V 



(2Xa7a + Xdc) 



(10) 



where H t is the heat capacity at constant stress, and "f a ,c 
are the Griineisen functions which describe the depen- 
dence of entropy on strain^ 

Using structural arguments, Moya et al. have asserted 
that both Xa and j a must be positive?^ This assertion in 
combination with Eq. (fTTJ|) suggests that the bulk ther- 
mal expansion coefficient a can be negative only if Xc 
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or 7 C is negative but not both. Early measurements 1 
yielded a negative value for Xc, which implied that the 
bulk thermal expansion coefficient of /3-eucryptite was 
negative because of the negative compressibility along the 
c-axis. Hortal et al^l put forth the idea that a negative 
Xc would explain the negative bulk thermal expansion in 
/3-eucryptite as an elastic effect associated with the in- 
terconnected Si and Al-tetrahedra as proposed by Gillery 
and Bushii£ 

Our results are, however, in contradiction to this point 
of view. The value we calculated for Xc using the elastic 
constants turns out to be positive. According to Eq. (fTU|) . 
a positive value for Xc implies that the Gruneisen func- 
tion 7 C must be negative in order to obtain a negative 
bulk expansion coefficient a. Indeed, recent phonon spec- 
tra calculations by Lichtenstein et al£& show that the 
Gruneisen parameters parallel to the c-axis for the modes 
around 400 cm -1 (bending of the Al-0 and Si-0 bonds) 
are large and negative, which leads to a negative value 
of 7c- Lichtenstein and coworkers attributed the nega- 
tive 7c to Li-position disordering and proposed an expla- 
nation for negative bulk CTE of /3-eucryptite similar to 
that of Schulzjiii Independently, Xu et al. used powder 
synchotron X-ray and neutron diffraction to show that 
cation disordering alters the structure of /3-eucryptite and 
significantly affects its thermal behaviour^ They have 
shown that Al/Si and Li disorder leads to a significant 
decrease in the lattice parameter c with only a moderate 
increase in a, leading to an overall volume contraction of 
~ l%.This behavior was explained^ as a combined effect 
of several interconnected phenomena including tetrahe- 
dral tilting, tetrahedra flattening, and shortening of the 
Si-0 and Al-0 bonds. 

Thus, our results are consistent with the findings of 
Lichtenstein et al£& and Xu et al£ leading us to con- 
clude that the negative coefficient of thermal expansion 
of /3-eucryptite is due to a negative value of 7 C associated 
with cation disordering, rather than to a negative Xc as 
proposed by Hortal et alJ^ 

V. CONCLUSION 

To summarize, we have computed the elastic stiffness 
constants of ordered /3-eucryptite containing 84 atoms 
per unit cell within the framework of generalized gradient 
approximation of DFT. The calculated elastic constants 



are in close agreement with the experimentally known 
values. The elastic constants were subsequently used 
to compute the linear compressibilities of /3-eucryptite 
parallel and perpendicular to the c-axis. Our calculated 
compressibility values agree well with those calculated 
from experimentally known elastic constants as reported 
by Haussiihl et al£2- The calculated values of compress- 
ibility are, however, in contradiction to those reported by 
Hortal et. al who measured the compressibilities Xc and 
Xa using a direct three-terminal method. Our calcula- 
tions show that the compressibility parallel to the c-axis 
is positive as opposed to the negative value obtained from 
the direct measurements ji£ 

Based on our calculations, we have also shown that the 
negative bulk thermal expansion of /3-eucryptite must be 
associated with a negative Griineissen function parallel to 
the c-axis rather than with a negative compressibility as 
proposed by Hortal et al. . The conclusion that the nega- 
tive bulk thermal expansion coefficient occurs because of 
a negative Griineissen function is consistent with the re- 
sults of Lichtenstein et al.fi who showed through the cal- 
culations of phonon density of states that the Griineissen 
function parallel to the c-axis is strongly negative due to 
the "bending" modes of the Si-0 and Al-0 bonds. Our 
results are also consistent with the neutron diffraction 
and X-ray synchrotron diffraction studies conducted by 
Xu et. al& 

The present study in conjunction with the results of 
Lichtenstein et al. and Xu et al. clearly indicates that the 
Xc must be positive and that the negative bulk thermal 
expansion is due to to cation disordering,- rather than to 
elastic effects^ 
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